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ABSTRACT 



Abstract: We present equations governing the way in which both the orbit and the 
intrinsic spins of stars in a close binary should evolve subject to a number of perturbing 
. forces, including the effect of a third body in a possibly inclined wider orbit. We 

I— i , illustrate the solutions in some binary-star and triple-star situations: tidal friction in a 



< 



wide but eccentric orbit of a radio pulsar about a B star (0045-7319), the Darwin and 
eccentricity instabilities in a more massive but shorter-period massive X-ray binary, and 
the interaction of tidal friction with Kozai cycles in a triple such as (5 Per, at an early 
stage in that star's life when all 3 components were ZAMS stars. We also attempt to 
• model in some detail the interesting triple system SS Lac, which stopped eclipsing in 

about 1950. We find that our model of SS Lac is quite constrained by the relatively 
good observational data of this system, and leads to a specific inclination (29°) of the 
outer orbit relative to the inner orbit at epoch zero (1912). Although the intrinsic spins 



o 



of the stars have little effect on the orbit, the converse is not true: the spin axes can 
vary their orientation relative to the close binary by up to 120° on a timescale of about 
Oh 1 a century. 
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Introduction 



We model the effect on a short-period binary-star orbit, and also on the spins of the two 
components, of the following perturbations: 

(a) a third body (treated as a point mass) in a longer-period orbit 

(b) the quadrupolar distortion of the stars due to their intrinsic spin 

(c) the further quadrupolar distortion due to their mutual gravity 

(d) tidal friction, in the equilibrium-tide approximation 

(e) General Relativity. 
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The third body's effect is treated only at the quadrupole level of approximation, although in 
principle it should not be difficult to go to a higher order if necessary. The stellar distortion terms 
are also treated only at the quadrupole level; it might be rather harder here to go to a higher 
order, but it is probably even less necessary. Each of the perturbing forces has been averaged over 
an approximately Keplerian orbit. We follow the analysis of Eggleton, Kiseleva & Hut (1998) - 
hereinafter EKH98 - for effects (b) - (d). The third-body perturbation comes from the same type 
of analysis (as does the familiar Schwarzschild- metric correction). 

We illustrate the model with some binary and triple examples: 

(i) the circularisation of an initially eccentric orbit of a neutron star around an obliquely-rotating 
massive normal B star, based on the SMC radio-pulsar 0045-7319 (Kaspi et al. 1994); we assume 
that the B star has spin inclined at a large angle (135°) to the orbit, and model the way in which 
the spin parallelises as well as pseudo-synchronises with the orbit 

(ii) the Darwin instability, i.e. the tendency for an orbit to de-synchronise if the spin angular 
momentum of a star is more than a third of the orbital angular momentum; and the eccentricity 
instability, in which rapid enough prograde rotation of the star causes the eccentricity to increase 

(iii) Kozai cycles, i.e. cyclic large-amplitude variation of the eccentricity of the inner binary binary 
due to a highly-inclined outer orbit, which in combination with tidal friction can make the inner 
orbit shrink even if it is initially too wide for tidal friction to be important. 

We base some of these illustrations on actual stellar systems, but the data on these systems is not 
(yet) sufficient to test the model rigorously. 

We then apply the model to the interesting 14.4 d binary SS Lac, which eclipsed for the first 
50 years of the 20th century, and stopped eclipsing later. Recently Torres & Stefanik (2000) - 
hereinafter TS00 - have demonstrated the existence of a third body in a ~ 700 d orbit, which if 
it is inclined at a suitable angle to the inner orbit could well account for the necessary change in 
the orientation of the orbit relative to the observer. We find that in order to accommodate the 
data given by TS00, we require the outer orbit to be inclined at about 29° (or 151°) to the inner 
orbit. We also require a specific value (37°) for the longitude of the outer orbit's axis relative to the 
inner orbital frame in 1912 (epoch zero). We find a modest variation with time in the eccentricity, 
and in the rate of change of inclination of the inner orbit to the line of sight. These cause us to 
revise slightly the masses found by TS00. Our model is fully constrained by the known data on this 
system, but it is not over-constrained, so that unfortunately we are not in a position to confirm 
that the model is correct. We can however make predictions for changes that should be capable of 
confirmation or refutation in about 20 years. 

In §2, we set out the equations governing the change in the orbit (and also in the rates of 
rotation of the two components), discussing the level of approximation that we use. In §3 we 
illustrate the behaviour of the equations in a number of straightforward cases, two without and one 
with a third body. In §4 we apply our model to SS Lac, and we conclude with a discussion in §5. 
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Some mathematical details are given in an Appendix. 



2. Equations for orbital evolution 

The evolution of the orbit under the influence of the perturbations listed in §1 is well expressed 
in terms of the following 5 vectors: e, the Laplace-Runge-Lenz vector, which points along the major 
axis in the direction of periastron and has magnitude e, the eccentricity; h, the orbital angular 
momentum vector, pointing perpendicular to the orbital plane and with magnitude h, the orbital 
angular momentum per unit reduced mass fx; q = h A e, which makes a right-handed orthogonal 
triad e, q, h with the previous two vectors, since e and h are always mutually perpendicular; and 
also the spin vectors Q . Q 2 of the two components. The vector q is along the latus rectum, the 
line through the focus parallel to the minor axis. It is also convenient to define unit vectors e, q, h, 
a right-handed orthogonal unit basis. This basis is not an inertial frame, of course, since it varies 
with time as the system evolves under the perturbations. But provided that the perturbations are 
sufficiently small that they do not affect the orbit by more than a small amount on timescales less 
than the period of the outer orbit, it is possible to estimate rather simply the rates of change of 
these 5 vectors in response to the 5 perturbative forces listed above. 

The equations (EKH98) governing the rates of change of e, h, Six, O2 can be written as follows: 

1 de — 

-— = (Z 1 + Z 2 + Z GR )q-(Y 1 + Y 2 )h-(V 1 + V 2 )e 
e dt 

-(l-e 2 ){5S eq e-(AS ee -S qq )q + S qh h} , (1) 

1 dh 

= Wi- 

+(l-e 2 )S qh e-(Ae 2 + l)S eh q + 5e 2 S eq h , (2) 



, = (Y 1 + Y 2 )e-(X 1 +X 2 )ci-{W 1 + W2)h 



h< ^t = MM-^ie + ^q + ^h} , (3) 

h lt = ^{-*2e + X 2 q- + W 2 h} . (4) 

An equation for q follows from differentiating the product q = hAe. The quantities Vi, Wj, Xj, Yi, Zi, 
one for each component of the inner pair, are given below. They arise from the quadrupolar distor- 
tion of the two components. The first two are dissipative terms, due to tidal friction, which tend to 
enforce orbital circularisation and synchronous rotation, at least in the absence of a third-body per- 
turbation. The next three are mainly non-dissipative perturbations due to quadrupole distortions, 
giving precession and apsidal motion; but X, Y do contain a small dissipative contribution which 
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tends to bring the stellar rotations into parallel with the orbit. The Z term, giving apsidal motion, 
contains a GR correction. The tensor S, with components S ee , S eq , ... in the e, q, h frame, is due 
only to the third body, and its components are also given somewhat further below - equations (15), 
(16). 



The dissipative terms V\, W\ are 



1 + f e 2 + f e 4 + <|e 6 llfi lh 1 



h 2 



8° 



Wi = 
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tpi 



(1 - e 2 ) 13 /2 



1 + f e 2 + f e 4 + A e 6 
(1 - e 2 ) 13 /2 



18w (1 - e 2 ) 5 



n lh 1 + 3e 2 + |e 4 



(1 



,2\5 



(5) 



(6) 



with similar expressions for *2. f2ih = Oi.h is the component of Q\ in the direction of h, i.e. 
parallel to the orbital axis, and oj is the mean angular velocity of the inner orbit, i.e. 2n/P. The 
tidal-friction timescale tp is estimated here, in terms of an inherent viscous timescale ty for each 
star, as 

1 9 Rf MM 2 1 

t\7i 



(7) 



tvi a 8 M 2 (1 - Qi) 2 

Mi, M2 are the masses of the two components of the inner binary, M their sum, Ri, R2 their radii, 
and a the semimajor axis. Qi is a coefficient measuring the quadrupolar deformability of the star; 
it is closely related to the apsidal motion constant (EKH98). For an n~3 polytrope, 



0.08M.R 2 



Q 



0.028 



(8) 



The intrinsic viscous timescale of the star, tyi, is not easily determined, but we use an estimate 
based on (a) the timescale on which the star would be turned over if most of the luminosity Li were 
carried by convection (Zahn Eggleton, 1977), and (b) a dimensionless factor 7 which comes from 
integrating over the star the rate-of-strain tensor (squared) of the time-dependent tidal velocity 
field: 

1 / T, \ 1/3 

71 7777^ , 7i~0.01 . (9) 



t V i ,x \3MiRfy 

The timescale tyi is of the order of years or decades. The quantity 7 - see Appendix - is determined 
by a model of the tidal amplitude as a function of radius through the star (EKH98), obtained by 
solving explicitly the velocity field required by the continuity equation if isobaric surfaces within the 
star are always to be equipotential surfaces - the basic assumption of the equilibrium-tide model. 

The contributions Xi, Yi, Z\ to the rotation of the axes due to rotational and tidal distortion 
of *1 (including the small contribution of tidal friction), are given by 

M 2 Ai ^ lh Q lc Q lq l + |e 2 + |e 4 
1 2/zu;a 5 (1 - e 2 ) 2 2^ F i (1 - e 2 ) 5 ' 1 j 
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M 2 A 1 n lh fl lq t fl lc l + fe 2 + |e 4 



Zi = 



2110)0? (1 - e 2 ) 2 2wt F i (1 - e 2 ) 5 
M 2 A X 



2fj,uja 5 



2Q 2 lh - Q 2 C - Q 2 q 15GM 2 1 + |e 2 + ±e 4 



+ 



2(1 -e 2 ) 2 a 3 (1-e 2 ) 5 



(12) 



Here /x is the reduced mass, and Oi e , Qi q are the components of Q\ in the directions of e, q. EKH98 
give the linear combinations Xi£l\ q — YiQ± e and Xi$li c + Y\Qi q , rather than X±,Yi directly: see 
the Appendix to this paper. The coefficient A\ is 

* = ^ ■ da) 

In all of equations (5) - (13), we interchange suffices 1 and 2 to find the corresponding term for the 
second component of the inner binary. 

The GR contribution to apsidal motion is 

ZGMu 

ZGR = ac 2 (l-e 2 ) • (14) 

The effect of a third body (mass M3) is included here only at the quadrupole level of approx- 
imation, following Kiseleva, Eggleton & Mikkola 1998 (hereinafter KEM98). At this level, the CG 
of the inner binary, and the third body, are unperturbed, and so the outer orbit is exactly constant. 
Like the inner orbit it can be described by a right-handed triad E, Q, H. Within the inner binary, 
there is a perturbative force which in the lowest approximation is linear in the vector separation: 
5fi oc Tijdj, where d is the separation of the inner pair. The tensor T depends on D, the separation 
of the outer pair (i.e. of the third body and the CG of the inner pair). Averaging T over an outer 
orbit, and then averaging the effect of 5f over the inner orbit, assuming that both orbits are only 
slowly varying on this timescale, we find that the tensor S of equations (1) and (2) is 

Sij = C (5ij — 2>HiHj) , 

C = M 3^out (lr>) 

4(M + M 3 ) W (l-e 2 ) 1 /2(i_ e 2 ut )3/2 • ^ ) 

|E| = e out is the eccentricity of the outer orbit, and tu ou t is the outer orbital frequency. Then the 
effect of the force Si on the vectors e, h are as indicated in equations (1) - (2), where on referring 
to the basis vectors e, q, h we have 

S ee = C(1-3H C H C ) , S eq = -3CH c H q , etc. (16) 

A somewhat surprising but welcome simplification is that after averaging T over D we have depen- 
dence only on H, and not on E, Q as well. Note that C is not a constant, because of its dependence 
on e. 
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Equations (1) - (4) are closed by the fact that a, uj, which appear in several places on the 
RHSs, are obtained in terms of e, h by way of the relations 

= r^v = ^ . (17) 



GM(1 - e 2 ) ' " V P J ° 3 

As e and h evolve away from their initial values, along with tli,Sl 2 , we have to compute various 
vector and scalar products: q = h A e, and Q e , ^q, H e , etc. 

The physical significance of V, W, X, Y, Z is perhaps most easily seen by splitting each of 
equations (1) and (2) into two pieces, one each for the moduli (e, h) and one each for the unit 
vectors e, h: 

1 A p 

- # = -Vi ~ V 2 - 5(1 - e 2 )S eq , (18) 

de 

— = {Z 1 + Z 2 + Z GR + (l-e 2 )(4S ee -S qq )}q 

-{Y 1+ Y 2 + (l-e 2 )S qh }h , (19) 
ITt = ~ W i- W 2 + ^ 2 S eq , (20) 

^ = {Y 1 + Y 2 + (l-e 2 )S qh }e 

-{X l +X 2 + {Ae 2 + l)5 ch }q . (21) 
Equations (19) and (21) can both be written as 

f = KAW , (22) 

K = {X Y + X 2 + X TB , Y 1 + Y 2 + Y tb , Z y + Z 2 + Z gr + Z TB ) . (23) 

Clearly K = (X, Y, Z) = Xe + Yq + Zh is the angular velocity of the e, q, h frame relative to an 
inertial frame. The terms with suffix TB arise from the third body, and can be readily identified 
with the corresponding 5-terms in equations (19) and (21). It is easy to see that q satisfies the 
same equation (22) as e, h. 

Equations (1) - (4) can be integrated numerically, using for example a four-stage Runge-Kutta 
procedure. However equation (1) as it stands has the slight problem, numerically, that in situations 
where e — > (usually as a result of tidal friction), e becomes undefined. This is easily solved by 
using instead equations (18) and (19). There is of course some redundancy in the e equation, but 
equations (18), (19) together are very well-behaved. There is not the same problem with equation 
(2), since h can hardly get to zero in realistic circumstances. Consequently, equations (2) - (4), 
with (18) and (19), i.e. 13 first-order ODEs in all, are quite readily integrated numerically as they 
stand. There are in fact two redundancies, since e. h = as well as e. e = 1. 

We find it convenient to use as our computational (and inertial) frame the initial (t = 0) 
orbital frame, say eo, q , ho- We need to be given a number of constant scalars and vectors, i.e. 
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Mj, Lj, Ij, for each of the inner pair of stars, and M3, H, e out for the third body and outer 
orbit. We also have to supply 13 initial values, for e, e, h, fii, L >- Some of the components are 
rather trivial, because by the above definition eo = (1,0,0) and ho = (0,0, 1) at t = — and of 
course q = (0, 1, 0). The quantities a, u at t = follow from e, h at t = by equations (17). 

The non-trivial initial vectors H, il > are all given directions in the obvious spherical-polar 
form, e.g. 

H = sin an cos /?h e + sin «h sin /?h q + cos oh h . (24) 

Thus an and flu are two of the three polar coordinates, the colatitude and longitude, of the vector 
H in the orbital frame. We view H as a vector starting at the focus, and intersecting a sphere which 
is centred on the focus and has North pole on the h axis. Longitude zero (on the equator) is on the 
e axis, i.e. the projection of periastron. The components of H are constant in the computational 
frame of eo, q , ho, but the components H Cl H q , which appear in equations (16) change with 
time because e, q, h change with time in the eo, q , ho frame, according to equation (22). 

The directions of Qi, S72 are given similarly, by pairs of angles aa 1 , (3^, etc. These angles 
have to be given initially, but of course the fl's, unlike H, vary both in the inertial frame and in 
the instantaneous orbital frame. 

In order to be able to determine the radial velocity curve and/or the eclipse light curve, and 
how they might change with time, we have to specify in addition the (constant) direction, J, from 
which the orbit is observed - constant (like H) in the eo, q~o, ho frame, but not of course in the 
e, q, h frame. We specify J by two angles in the same way as H, and call them aj, (3j. The angle 
a j is just the usual inclination of the orbit to the line of sight. The angle j3j is almost the same as 
the 'longitude of periastron'. The latter quantity is usually called oj, but we call it lu\ p as we have 
already used u> for the orbital frequency. The relation between (3j and u\ p is 

/3j+w lp = 270° . (25) 

We think of aj, (3j firstly as given initial conditions; but at later times they can be evaluated from 

cosaj = J.h , tan/3j = £^ . (26) 

J. e 

J is a constant in space, but variable in the orbital frame e, q, h since these unit vectors vary 
with time. The same formulae (26), mutatis mutandis, give the corresponding a (colatitude) and 
(3 (longitude) for the other vectors H, fli, and O2 at later times. 

If we are fortunate enough to have very accurate observations over sufficiently long stretches of 
time, we may be able to measure some rates of change such as P, e, cij, /3j. Equation (18) already 
gives e, and along with equations (17) and (20) this gives P (or u or d): 
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The third-body terms cancel, because they are conservative and do no work around a Keplerian 
orbit at our level of approximation: only tidal- friction terms affect P (or a). We can differentiate 
equations (26) w.r.t. time, keeping J constant and using equation (22) for any of e, q, h: 

J. K A h ■ J.h J.K-K.h 

| J A h| |JAh| 2 

K being the rotation rate of the frame as given in equation (23). The second relation involved some 
elementary vector manipulations, but the first came directly from equations (22) and (26), using 
sinaj = | J A h|. 

We wish to emphasise the following point, which we believe is treated incorrectly in much of 
the literature which we have read. The rate of rotation of the line of apses, i0\ p = — /3j, is usually 
attributed to Z, the h component of K, to the extent that Z is normally referred to as 'apsidal 
motion'. But it is easy to see that (3j in equation (28) can be non-zero on account of the e, q 
components X, Y of K as well, as will happen with a massive rotating star whose spin axis is 
highly inclined to the orbital axis and precessing about it. We can see that if X, Y = 0, so that 
K = Zh, then uj\ p = — (3j = Z as expected. But if X, Y are not zero (precession) they contribute 
to /3j, even in the case that Z = 0. Note that this effect does not depend in any way on the details 
of our model: it only depends on the fact that the orbital frame e, q, h has some general angular 
velocity K = Xe + Yq + Zh, while the system is viewed from a fixed direction J with variable 
colatitude aj and longitude f)j in the e, q, h frame. 

The effect of precession on f3j is mainly to swing the line of apses back and forwards (libration) , 
rather than to advance it monotonically (circulation), as does Z. But for the case where (a) is 
not parallel to h, (b) the tidal friction terms in equations (10) and (11) are negligible - which they 
usually are - and (c) the ^-dependent term in equation (12) dominates over the remaining term - 
which is usually the case for rapidly-rotating components - the librating and circulating effects are 
comparable. 

Our model can be used to provide times (or phases) of eclipses, when the inclination is sufficient 
for eclipses to occur. We use the simplest approximation, that the stars are spherical. For the 
beginning and end of an eclipse we have to satisfy the equation 

| J A d| = Ri + R 2 ■ (29) 

Here d, the vectorial separation of the two stars, has components in the e, q, h frame given by the 
usual formula 

d = (ecosfl + qsinfl) , (30) 

1 + e cos 6 

with / = a(l — e 2 ) = h 2 /GM being the semi-latus-rectum. Equations (29) and (30) give a quartic 
equation for cos 6, which we solve analytically. The coefficients of the quartic are determined by 
J.e and J.q, which vary with time as the basis set moves. Having determined the four, two or zero 
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real roots that lie in the range [-1,1], we can determine the phase <E> (i.e. time, divided by period) 
of ingress and egress from the usual formulae 

e ~t~ cos 9 

2tt$> = ib — sin ib , cosip = . (31) 

l + ecos# v ' 

This is phase measured from periastron (9 = 0). We can also determine the phases of other 
significant points on the orbit. The point where the projection on to the orbital plane of the the 
line-of-sight vector J intersects the orbit (conjunction) is given by 6 = 0\ say, where 

(J — J. hh) A d = 0, i.e. tanfli = =5 . (32) 

J.e 

The point where the radial component of velocity (relative to the CG) vanishes is given by 9 = #2 
say, where 

J.d = 0, d = j(-e sine + {e + cos 0}q) , (33) 

so that 

sin(0 2 -0i) = e sin 0i . (34) 

Thus the values of 9 at both these points are also functions of J.e and J.q. There are two values 
of both 9\ and 9<i in the range — 360°; in §4 we take 9\ such that *1 is in front, and 02 such that 
*1 is behind. For triple systems we ignore the small effect due to the variable motion of the CG 
of the inner pair. The points where the radial velocity is a maximum (or minimum) are given by 
J.d = (since d || d in the unperturbed orbit), and hence by #3 = 9\ ± 180°. 

The model we present here has, we believe, the merit of considerable simplicity, both concep- 
tually and numerically. We emphasise here the approximations on which it is based: 

(i) Only the quadrupolar component of the distortion of each star is modelled. This assumption 
may be fairly good in systems which are only mildly eccentric, but can be expected to be less valid 
in systems of high eccentricity. For the distortion due to rotation, it is assumed that the stars are 
in solid-body rotation. 

(ii) The components are assumed to adjust instantaneously to fill an equipotential of the joint 
gravitational-centrifugal potential. This leads to a specific tidal velocity field within each star, whose 
shear, combined with a viscosity assumed to be due to convectively-driven turbulence, determines 
the force of tidal friction. Although some analyses have argued that the effect of convection in 
a star with a radiative envelope and a convective core is small, we follow the analysis of EKH98 
which shows that tidal dissipation within a convective core is not small. 

(iii) The effect of the third body is only modelled at the quadrupole level of approximation. This 
is sufficient to demonstrate such phenomena as Kozai cycles (Kozai 1962), where the third body, if 
placed in an orbit highly inclined to that of the first two, causes large fluctuations in eccentricity on 
a long timescale. The approximation is not very good for a 'parallel orbit, since all the off-diagonal 
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components of the tensor S vanish so that the only effect in equations (18) - (21) is the apsidal- 
motion term. It therefore does not allow us to model the fluctuations in eccentricity that the third 
body produces in the inner pair; actually these are quite small, according to an N-body integration, 
but they can be significant on a long timescale when combined with tidal friction (KEM98). 

(iv) The same level of approximation means that there is no additional force, or couple, on the 
outer binary. Consequently angular momentum is not conserved: the inner binary can gain or 
lose angular momentum, but not the outer. Moderate accuracy relies on the fact that the angular 
momentum of the outer binary is large compared with the inner, so that a substantial amount 
lost by the inner counts as only a small perturbation to the outer. However angular momentum 
increases with only the cube root of the period, at given masses, so that a period ratio of 50 (an 
unusually small value, but appropriate to SS Lac) means an angular momentum ratio that is not 
large. 

3. Some illustrative examples 

Our model for perturbed orbits is original to the extent that (a) it has a specific formulation 
for the parallelisation of stellar spin which is initially oblique to, or even anti-parallel to, orbital 
angular momentum, and (b) it includes the effect of a third body along with the other perturbations. 
Soderhjelm (1975) gave a formulation of the third-body effect, but without the other effects. As 
applied to systems which are binary rather than triple, and where the stellar spins are (at least 
by hypothesis) parallel or nearly parallel to the orbit, our model does not differ from lowest-order 
standard analyses. We confirm the following standard results: 

(i) on a timescale of ~ixF I£l/fJ,h, the spin becomes parallel to the orbit and 'pseudo-synchronous', 
i.e. it reaches the value where the viscous couple W in equation (6) is close to zero (Hut 1981). 
The couple is an average around a Keplerian orbit, and it vanishes when the larger but short-lived 
couple near periastron is balanced by the weaker but longer-lived (and opposed) couple at apastron. 
Equating W to zero gives the pseudo-synchronous value of J\ as u> multiplied by a function of e. 

(ii) On the longer timescale £tf the orbit is circularised. 

However, both these statements have to be qualified by the condition that the spin angular momenta 
of the stars have to be suitably small when compared with the orbital angular momentum; otherwise 
the binary can become either desynchronised or decircularised. 

(iii) For triples, ignoring the effects of quadrupolar distortion, tidal friction and GR, we obtain 
equations which allow the eccentricity and the mutual inclination of the inner orbit to fluctuate 
periodically between limits (Kozai cycles: Kozai 1962, Mazeh & Shaham 1979; KEM98). If we start 
with e = and sinan > \/2/5 («h~39°), these cycles can have large amplitude. The maximum 
eccentricity reached is 

e max = f sin2 «H - § ~ 1 - §<5a H 2 if <5a H = ^ - a H . (35) 
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We see that e max can approach very close to unity if «h is only moderately close to 90°. The 
timescale of these cycles is of order 1/C - equation (15) - and so of order P^ nt /P, apart from a 
mass-ratio-dependent factor which is only significant if the third body is much less massive than 
the other two. 

Commonly, among observed binaries, either the orbital period is sufficiently short that tidal 
friction has already circularised it, or sufficiently long that tidal friction is insignificant on a nuclear 
timescale. This is because of the high power of R\/a in equation (7). There is only a fairly narrow 
range of periods, perhaps 4 — 5 d (but depending on mass and age) where one might hope to find 
binaries in which parallelisation is still taking place. However, there exists the interesting SMC 
radio-pulsar binary 0045-7319 (Kaspi et al. 1994) in which it is conjectured that, as a result of an 
asymmetric supernova explosion (SNEX), the neutron star is on an inclined orbit, relative to the 
spin of the normal B1V component. In our first example below, we endeavour to model the process 
of parallelisation etc. of the B star in this system. In our second example, on the supposition that 
an SNEX may typically put a neutron star (NS) into a non-synchronous orbit, we also consider a 
model of a massive star with an NS companion, and various degrees of asynchronism. Such models 
can experience desynchronisation and/or decircularisation. 

There is only a rather small number of known triples where the inclination of the outer to 
the inner orbit is directly measured, and even fewer in which it is clearly established that this 
inclination is large enough to cause Kozai cycles. In fact, the system (3 Per (Lestrade et al. 1993) 
is the only example we know. We therefore consider how the orbit of this system might have been 
modified at an early stage in its life, when it was a detached near-ZAMS system. 

In binary orbits that are eccentric, but already (at least by hypothesis) parallelised and pseudo- 
synchronised, the only part of our model to be testable is the effect of tidal distortion and GR on 
apsidal motion. In this respect our model is no different from the classical model: Claret & Gimenez 
(1993) have discussed apsidal motion, comparing observed values with those expected theoretically 
from the combination of quadrupole distortion and GR. For many systems there is reasonably good 
agreement. For some systems however there is disagreement, even strong disagreement, for example 
DI Her (Guinan et al. 1994) and V541 Cyg (Lacy 1998). We suspect that the discrepancies here 
may be due to the presence of a third body, so far undetected; although another possible reason 
for aberrant apsidal motion is that the stars are rotating obliquely to the orbit. The latter leads 
to smaller apsidal motion than expected for parallel synchronism, and the apsidal motion can even 
have the opposite sign - equation (12) - if ^ Qi/y/3. However a third body can also contribute 
apsidal motion of either sign. 

3.1. (i) Parallelisation, synchronisation and circularisation in a wide eccentric orbit. 

We consider the effect of the perturbative forces within a binary roughly based on the radio- 
pulsar binary 0045-7319 in the SMC (Kaspi et al. 1994, Bell et al. 1995). In this binary (no third 
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body is detected, or suspected) the pulsar's orbit is quite wide (P = 51.17 d), and highly eccentric 
(e = 0.808). The directly-measured longitude of periastron gives f3j = 270° — uj\ p = 154.76°, from 
equation (25). The companion B1V star appears to be unusually inactive: it is not a Be star, 
apparently has negligible wind, and the pulsar is not accreting significantly, at least not enough 
to be an XR source. These unusual circumstances (for massive neutron-star binaries) mean that 
the radio orbit is unusually well-defined, so that even the slight change of orbital parameters, due 
presumably to tidal friction, apsidal motion and precession, are measurable. 

We refer to the NS component as *1, because it is descended from what was presumably the 
originally more massive component, and the B1V star as *2. This choice determines which suffix 
belongs to which star. 

Although the mass-function of the pulsar is accurately known, there is only a very tentative 
radial-velocity curve for the B star. Bell et al. (1995), assuming Mi = 1.4 Mq, suggest the follow- 
ing parameters: M2~8.8M , i?2~6.4i? , L2~1.2 x 10 4 L Q , with substantial uncertainties. The 
consequential inclination of the orbit to the line of sight is aj ~44° (or 136°). Bell et al. also esti- 
mate a projected rotational velocity for the B1V star of V TOt sin % = i?2|^2 A J| ~ 113km/s, which 
suggests a rotational period for the star of less than three days, but depending on the unknown 
orientation of the stellar spin relative to the observer. The spin rate though not clearly known 
is marginally consistent with the possibility that the B star is in pseudo-synchronism (Hut 1981). 
For e~0.8 pseudo-synchronism requires f^h ~ 12.5a; or i-*~4d. It is very likely however that the 
NS was put into its current highly eccentric orbit by a supernova 'kick', which also makes it likely 
that the stellar spin is inclined, perhaps quite highly inclined, to the orbit. It is in fact easier to 
account for the rate of orbital period change if the spin is retrograde, since this tends to maximise 
the contributions of W2 and V2 in equation (26). 

Kaspi et al. (1996) further determined various rates of change: P/P 2.2 x 10~ 6 /yr, dj = 

2.1 x 10~ 4 rad/yr, = -4.5 x 10~ 4 rad/yr. The sign of dj will be different if we adopt aj ~ 136° 
instead of 44°. The accuracy of these quantities, and V vot sini, is of the order of 3 — 10%. 

Our model is slightly over-constrained by the current observational data, supposing that we 
take literally the estimate (9) for the viscous timescale. We adopt the values of M2, R2,L2, Mi, P, e, aj, (3j 
mentioned above. There is no evidence for a third body, and so we take M3 = 0. We ignore all 
parameters relating to the neutron star except its mass, since its spin angular momentum will be 
too small to influence the system. This only leaves the three components of Q2 to be assigned at 
t = 0, and there are four constraints to be satisfied: P/P, V vo tsm.i, dj and j3j should all have the 
values listed above. 

Fig 1 is a short evolutionary run starting from Q/oj = 20, an 2 = 135°, /3q 2 = 0°. On such a 
short timescale only /3q 2 changes significantly, by precession. The plotted quantities are the ratios 
of the computed to observed values for the four quantities listed. It can be seen that at 195 yr all 
four ratios are fairly close to unity, at which point /3q 2 = 111°. We therefore adopt this as a new 
starting value. 
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FIG 1 

Since the starting values used above were a shot in the dark, we can expect to get better 
agreement by some procedure such as least-squares. However, the answers will be very strongly 
dependent on (a) the radius R2, which enters to the eighth power in equation (7), and (b) the 
estimate (9) for 7, which must be very uncertain. It might be more realistic to treat 7 as an 
unknown, in which case probably an exact solution (and possibly several, because of the non- 
linearity of the equations) can be found; but it will still be very dependent on R2, which cannot 
be accurately known. Thus we feel it is premature to attempt a definitive solution, but we feel 
encouraged by the fact that the model is not obviously wrong. 

The evolution of the eccentricity, of the component of spin parallel to the orbit (relative to the 
total spin), and of the orbital period (relative to initial period) is shown in Fig 2a, for a timespan 
of about 3Myr into the future. We started with the parameters listed above (but /3q 2 = 111°). 
The evolution of the B star in this interval has not been allowed for: the main-sequence lifetime 
of an 8.8 M & star is expected to be about 33 Myr. The perpendicular spin goes through zero at 
about 0.5 Myr, and the spin is almost completely parallel by 1.7 Myr. Circularisation takes a good 
deal longer, and is only half-complete by 3 Myr - but it will start to be strongly accelerated by the 
neglected evolutionary expansion, at this stage. Currently e/e~ — 2.5 x 10" 7 /yr, on our model. 
This is comfortably below the upper limit found by Kaspi et al. (1996) of 7 x lCT 6 /yr. 

FIG 2 

Fig 2b shows the two components of spin in the orbital plane, £22 e and ^2-q> plotted against 
each other. To make the figure clearer the viscous evolution was speeded up by a factor of 10 8 / 3 ; this 
makes for a much less tight spiral pattern. Evolution starts slightly left-of-centre at the top edge. 
The rotation axis precesses counter-clockwise about 1.25 times, until the vertical component of O2 
(02-h, Fig 2a) changes sign, and then precesses clockwise while the two horizontal spin-components 
diminish towards zero. Had we kept to the more realistic viscous timescale of Fig 2a, there would 
have beeen about 550 revolutions of the axis before it reversed direction. 

Fig 2c shows the evolution of four timescales, also using the speeded-up model of Fig 2b. The 
timescales are all given as logs, and in years. The timescale for period change \P/P\ (plusses) starts 
at ~ 10 31 yr, which would be roughly the required value of ~5 x 10 5 yr if we did not speed up the 
viscous evolution by 10 8//3 . The eccentricity timescale |e/e| (asterisks) is about 6 times longer to 
start with, but is more nearly constant. The two other timescales shown are both related to apsidal 
motion: \/Z (circles), and (equation 27; crosses). /3j is the actual apsidal motion, which 

however is influenced by the precessional terms X, Y as well as by the usual 'apsidal motion' term 
Z. Prior to about 1000 yr, when the vertical spin passes through zero in the speeded-up model, the 
line of apses turns at a highly variable rate; probably the axis was librating rather than circulating, 
in the fairly recent past. Once the B star is no longer counter-rotating the line of apses circulates 
more uniformly, but with an oscillating component which diminishes as the spin becomes more 
parallelised. 
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Figure 2d shows two more timescales (also logged): the precessional timescale, 1/(X 2 +Y 2 )~ 1 / 2 
(plusses), and the the timescale of the change of inclination of the orbit to the line of sight, l/|dj| 
(asterisks). The precessional rate goes through zero once, causing the cusp at ~1000yr. The 
orbital direction oscillates about zero, causing many cusps in the log modulus of its derivative. 

The origin of the present system presents some puzzles, and has been the subject of recent 
controversy (van den Heuvel & van Paradijs 1997; Iben & Tutukov 1998; hereinafter HP, IT). HP 
favour a history that involved Roche-Lobe overflow (RLOF) followed by a supernova explosion 
(SNEX) with an asymmetric kick, and IT a history that involved a common-envelope (CE) phase 
followed by an SNEX without a kick. We believe that neither history is satisfactory, and propose a 
scenario which is somewhat similar to IT in its earlier phase (but requiring a less massive progenitor 
to the NS), and rather like HP in its later phase, requiring an SNEX kick. 

We would normally expect that the system, having started with two massive MS stars, would 
have evolved through RLOF, so that *1 (the originally more massive component) would have 
become a helium star, perhaps with a modest H-rich envelope, before heading on to C-burning and 
so fairly quickly to an SNEX (HP). However two things argue against this: 

(i) If *1 was originally over ~ 10 M@, enough to leave a post-RLOF remnant capable of an SNEX, 
then RLOF should have made *2 considerably more massive than it is now (even though its mass 
is by no means certain). In such RLOF we normally expect *2 to become more massive than the 
original mass of *1. 

(ii) We would expect that as a second result of the RLOF *2 would be a very rapid rotator, a Be 
star more-or-less, instead of the rather slowly rotating and unusually inactive star observed. 

A possible answer to both these points is that 

(a) the initial *1 was only moderately more massive than *2 now, say k, 12 M Q 

(b) the initial *2 was little different from the *2 now seen (the Bl star) 

(c) the binary was fairly wide initially, say P k, 50 d 

(d) *1 evolved to a point where its outer layers, helped by the disturbing effect of the binary 
companion, became unstable and blew away, firstly as a P Cyg star and then as a Wolf-Rayet star, 
perhaps without *1 ever reaching a radius as large as its Roche-lobe radius. 

If *1 did reach RLOF, this might have been more like a CE event, with much of the envelope 
disappearing to infinity rather rapidly, and with only moderate, or perhaps even negligible, orbital 
shrinkage. But we would rather categorise the process as 'binary-enhaced stellar wind' (BESW), 
which may have altogether prevented *1 from ever reaching its Roche lobe. 

The WR binary 7 2 Vel, with P = 78.5 d and e = 0.33 (Schmutz et al. 1997) might be 
of the same character as the possible immediate precursor to the 0045-7319 binary, since the high 
eccentricity argues against there having been an episode of Roche-lobe overflow. The *2 of j 2 Vel 
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is an 08III star of 20 — 30 M Q , substantially more massive than we require. Consequently *1 would 
also have been substantially more massive originally, perhaps by about the same factor. 

IT's model was somewhat similar to ours, except that they argued for a more massive initial 
*1, ~28M . This was required because of their desire to produce the configuration of 0045-7319 
without an SN kick. They argued for a CE event which reduced the period from an initial value of 
27 — 76 d to a value of ~3.2d. They postulate that the obliquity of the spin to the orbit, strongly 
suggested by the measured dj of Kaspi et al. (1996), is simply left over from a primordial obliquity, 
and survived any possible tidal friction during the helium-star phase, when in their model the 
orbital period was 3.2 d. A difficulty with this is that with *1 initially so much more massive than 
*2, *2 should be rather little evolved, and should be substantially smaller than the value of 6.4 R@ 
suggested by Bell et al. (1995). Our model supposes a much less massive *1, and so allows *2 to 
be more substantially evolved. Our model does not predict, nor need to predict, the orbital period 
during the helium-burning phase; we accept the probability of an asymmetric kick, which could 
in principle lead to the present period if the intermediate period was anywhere in the range of 
~3-50d. 

Our model of the current orbital evolution might give an upper limit to the age of the system 
(since the SNEX), by integrating backwards from present conditions. This is not a very safe process, 
numerically, in a dissipative system, but we made an estimate of the accuracy by integrating 
forwards again. It appears that in fact the evolution decelerates going backwards, as is hinted 
at by the behaviour of Q,2h/^2 in Fig 2a. We integrated back ~9 x 10 5 yr, reaching P = 200 d, 
e = 0.922 and ajj 2 = 149°; on integrating forwards again we recovered P, e, an 2 and f^/u; to 5 
significant figures, while /3q 2 was in error by about ~80° after several thousand rotations of the SI2 
axis. The spin period ~ 9 x 10 5 yr ago is predicted to have been 1.6 d. Although it is marginal, this 
may be consistent with the B1V component's still being reasonably inactive, so that the model is 
still applicable. Thus it is possible that the system may be as much as ~ 10 6 yr old in its present 
form. The required orbit so long ago might seem improbably long and eccentric; but one might 
reasonably think the present orbit improbably long and eccentric if it had been hypothesised rather 
than measured. 

The equilibrium-tide model of tidal friction has often been considered inadequate for systems 
like 0045-7319. There appear to be two main reasons, one of which we largely accept and the other 
we reject. In order, they are 

(a) near-equilibrium is not very likely to be established in a highly eccentric orbit; it is more 
reasonable in a nearly circular orbit (like the Earth-Moon system) 

(b) Although turbulent convection may be a good source of friction in stars with deep convective 
envelopes, massive stars are only convective in their cores where the amplitude of the tide is 
considered to be too small to be significant. Radiative damping in the outer layers might contribute, 
but this is orders of magnitude smaller. 

We believe that (b) is largely based on a highly inexact estimate of the equilibrium-tide velocity 
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field. 



If the principle is accepted that surfaces of constant density (and pressure) are always closely 
equal to equipotential surfaces (the basic assumption of the equilibrium-tide model) then presum- 
ably the velocity field is determinate, and comes basically from conservation. Alexander (1973) and 
Zahn (1977, 1978) made crude estimates, based on the motion being assumed either incompressible 
or irrotational, and concluded that the amplitude of the tide (whose square is proportional to the 
rate of dissipation) goes to zero like r 4 , approaching the centre. If the convective core were, say, 
one third of the stellar radius, then the dissipation would be down by ~ 10~ 4 relative to a star 
with a largely convective envelope. However, EKH98 determined - their equations (100) to (112) 
- an expression for the tidal velocity field, and its rate of viscous dissipation, which is exact, to the 
extent that (a) the equilibrium-tide model is exact, and (b) dissipation is primarily by the effective 
viscosity of turbulent eddies. The velocity field is neither irrotational nor incompressible, nor does 
it diminish to zero like r 4 . Rather, the tidal amplitude diminishes from its surface value by less 
than a factor of 10 for typical MS models. Thus the effect of dissipation in the convective core is 
by no means negligible: it may be down by ~ 10~ 2 only. This is the basis for our estimate of 7 in 
equation (9). In the Appendix, we briefly summarise the analysis of EKH98 regarding the factor 
7- 

Witte & Savonije (1999) computed the spectrum, and damping rates, of normal modes that 
can be expected to be excited in a 10 M & star as a result of perturbation by an NS companion 
with the orbital parameters of 0045-7319. They obtained a braking timescale that was usually in 
the range 10 55 — 10 6 ' 5 yr. The timescale changes rapidly by factors of ~ 10, both up and down, on 
timescales of only 10 4 yr or less. There are occasional excursions to values of braking timescale as 
low as 10 3 yr, which result from modes resonating with harmonics of the orbital frequency. There 
are also occasional episodes of orbital spin-up rather than spin-down. Such a detailed model may 
well be demanded by the physics, but inevitably means that the interior structure of the star will 
have to be very precisely known: a good deal more precisely than information which is currently 
available. We hope that our estimate, equation (7), can serve as a crude average over a range of 
time of more detailed values that can only be computed if the structure and rotation of the star 
are known to considerable accuracy. 



Inherent in equations (1) - (4) are at least two kinds of instability. Firstly, there is the Darwin 
instability. Consider the case of a binary (i.e. no third body) where the spin of one massive 
component (*2) is parallel to the orbit, and the companion (*1) is a point-mass neutron star, as in 



3.2. (ii) The Darwin and eccentricity instabilities 



the previous example. Equation (25) for to can be united with equation (4) for 0,2 to give 
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' lh fa(e)-3f b (e) 



^2 



/c(e)-3/ d (c) 



(36) 



where the functions f a ,---,fd are all functions of e that can be evaluated from equations (5) and 
(6). All these functions tend to unity as e — > 0. It can be seen that as long as 

^ - 3/d(6) - 3 if e ~ , (37) 



m 2 f c (e) 

then 0,2 ^ as time increases. But if inequality (37) is violated, f^/w will diverge as time 
increases. For a value of e which is not small, there still is a critical condition but it is e-dependent. 
This well-known instability reqires that the spin angular momentum must be greater than a 
third of the orbital angular momentum fih (for e = 0). 

The eccentricity instability is seen in equations (18) and (5). Specialising once again to the 
situation where one star is a point mass, and e = 0, we see that if 

18 

n 2h > — lo , (38) 

then the eccentricity starts growing exponentially. If e > to start start with, there is still the 
possibility of e growing, although the criterion is now e-dependent. In other words, if the star is 
rotating fast enough, it gives up its angular momentum in spurts sufficiently concentrated towards 
periastron that the companion star is flung into a wider and wider orbit - but with periastron not 
much changed because that is where the largely tangential impulse peaks. 

FIG 3 

Although both instabilities give exponential growth the result can sometimes be surprisingly 
self-limiting. Fig 3a shows the evolution of a system whose initial configuration was unstable to both 
processes. We took a very massive star (40 M ) evolved substantially across the MS (to 20 R & ), 
put it in a 6d, e = 0.1, orbit with a neutron star of 1.4 M Q , and started it in parallel rotation 
at twice the orbital rate. We used the default values (8) of moment of inertia and quadrupolar 
distortion. Both the eccentricity and the degree of non-corotation (measured by Q,/uj) began to 
grow. For stars of comparable mass it is difficult to violate criterion (37), but if one star is much 
more massive than the other, it can also be large enough, without quite filling its Roche lobe, to be 
Darwin-unstable. However although the star spins up relative to the orbit, the orbit gains angular 
momentum, and so loses angular velocity, fast enough for the D-stable criterion (37) to become 
satisfied later. After between 10 6 and 10 yrs, the orbit first becomes D-stable and later E-stable, 
and tends to both synchronism and circularity with a period of ~30d. However, as before we have 
ignored nuclear evolution in the massive component, which would no doubt fill its Roche lobe in 
little more than 10 6 yr. 

Fig 3b is the same system except that the stellar spin rate was started at 70% of corotation, 
rather than twice. This is substantially E-stable and very slightly D-stable to start with, but as 
e decreases towards zero, and U/u increases (though only very slightly) towards unity, at about 
5000 yrs the system crosses the D-unstable margin. Although both Q and u) are going up, trying to 
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approach synchronism, while IQ obviously goes up \ih goes down because the orbit shrinks. Hence 
the D-stable criterion (37) crosses into instability and the system begins to move rapidly away from 
corotation. The orbit continues to circularise, but desynchronises and shrinks rapidly towards a 
collision at ~ 19000 yr. 



3.3. (iii) Kozai cycles with tidal friction. 

We now consider a problem which has a third body as well as quadrupole distortion and 
tidal friction. When the outer binary is sufficiently inclined to the inner binary it is possible for 
the eccentricity of the inner binary to fluctuate slowly by a large amount (Kozai cycles). The 
amplitude of the eccentricity fluctuation depends only on the inclination, and not on the outer 
period, or eccentricity, or third-body mass; the period of the fluctuation is of order 2ir\fl — e 2 /C 
- equation (15). Even if the inner binary, when it is nearly circular, is too wide for tidal friction to 
play a role, the increase in eccentricity may make tidal friction important at some point in the Kozai 
cycle. Recall that a, like to and P, is unaffected by the third body at our level of approximation, 
as shown by equation (27), so that as e increases the periastron separation decreases. We illustrate 
this with the well-known semidetached binary Algol {(3 Per), which has a third body(~ 1.7 M©) in 
a 679 d orbit inclined at 100° to the semidetached pair's orbit (Lestrade et al. 1993). 

In its present configuration, the inner pair is not subject to Kozai cycles, because the per- 
turbation due to the quadrupole distortion of the lobe-filling component is much larger than the 
perturbation due to the third body. However, at an early stage in its life (5 Per must have been a de- 
tached binary of two near-ZAMS stars, with radii and therefore quadrupole moments substantially 
smaller than at present. 

If we believe that (3 Per has evolved without mass loss (ML) or angular momentum loss (AML), 
i.e. conservatively, we would be able to infer the period at any mass ratio, from 



(1+qf = Mi 

q 3 ' q ~ M 2 



P oc v±Z*L_ , q = £± . (39) 



Taking an illustrative q$ = 1.25, the present period P = 2.87 d and q = 0.216 imply that Pq ~0.6d. 
However, although we accept provisionally that ML may have been negligible, there is direct and 
indirect evidence that cool Algols experience AML, presumably by magnetic braking in a stellar 
wind (Refsdal, Roth & Weigert 1974, Eggleton 2000a). For given masses, the period goes like h 3 , 
and so if the system lost 50% of its angular momentum, it must have started with Po ~4.8d. 

What we show in this subsection is that the initial period, had it been longer than ~ 3 d, would 
have shrunk, by a combination of Kozai cycles and tidal friction, to a value under 3 d in a fairly 
short interval of time ( ^ 10 7 yr). Consequently we have an upper limit to the amount of AML that 
could have taken place subsequently, once *1 became a cool subgiant subject to magnetic braking 
(Eggleton 2000b): about 40% of the initial angular momentum. 
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Fig 4 shows the evolution of a 'proto- Algol' system with an initial period of 5 d, in (a) the short 
term, (b) the medium term, and (c) the long term. The initial Kozai cycles reach up to e = 0.67 
(starting somewhat arbitrarily at e = 0.1). This value is well short of the maximum that would be 
reached (e = 0.985) if quadrupolar distortion was negligible, but is nevertheless quite large. Tidal 
friction near periastron at the peak of the Kozai cycles reduces the range of variation of e, though 
somewhat unexpectedly by increasing the minimum eccentricity even more than by reducing the 
maximum. By about 10 6 yr the eccentricity fluctuates between 0.47 and 0.53, and both the range 
and the mean reduce until by 10 7 yr the orbit is circularised at P~2.1d. 

FIG 4 

A point to note is that the inclination an of the inner orbit to the outer orbit changes somewhat 
during the process. We started from 97.5°, in order to end up with the currently observed value of 
100°. For longer initial periods the change is larger, which probably means that the period was not 
in practice much larger than ~ 10 d before the Kozai cycling and tidal friction reduced the period 
to ~2d. 

It is not clear how triple systems, and especially such close triple systems, formed in the first 
place, but a possible mechanism, arguably the least unlikely, is that, fairly early on in the star- 
forming process when the stellar density was higher than it is now, two primordial binaries had 
a near-collision, with one component of one binary captured by the other binary, and the other 
component ejected. In this scenario, angles near 90° are much more likely than those near 0°. 

Table 1 shows how the the period P en d at the end of the shrinkage process depends on the 
period Pq at the beginning, for our specific proto- Algol system. It also shows the time taken in the 
shrinkage and circularisation process, which is always small compared with the expected nuclear 
lifetime of the system ( ~ 1 Gyr) , and gives the starting value of mutual orbital inclination «h that 
will end up as the current value of 100°. Assuming that this inclination is distributed randomly, in 
a capture process, the range 80 — 100° has probability ~ 17%, and the range 86 — 94° about 7%. 

TABLE 1 

The combination of Kozai cycles plus tidal friction should mean that there is a shortage of 
triple systems with (a) fairly high inclination of one orbit to the other, and (b) inner periods above 
perhaps 3 - 4d. This will be difficult to confirm, because it is very difficult to determine the 
inclination of one orbit to another. SS Lac (below) is a triple in which we infer an~29°, which 
is not enough to give significant Kozai cycling; thus the inner period of 14 d does not conflict 
with our conclusion. If inclinations are indeed distributed at random, then ~ 50% of triples have 
60° & a j & 120°, and a quite substantial deficit of systems with inner periods longer than say 3 — 4d 
can be expected. 
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4. An application to SS Lac 

SS Lac is a binary which eclipsed before about 1950, but not subsequently. A likely explanation 
was the presence of a third body, in a non-coplanar orbit, and this was confirmed by Torres & 
Stefanik (2000) - hereinafter TS00 - who found long-period orbital motion in the CG of the short- 
period pair. By coincidence, the longer period in SS Lac is exactly the same as that in Algol (679 d). 
Following TS00, we refer to the 3 components as Aa (*1), Ab (*2) and B (*3), and the two binaries 
as A and AB. TS00 also re-analysed historic light curves of the period 1890 - 1930, obtaining a 
mean light curve assigned to epoch 1912. Their spectroscopic data refers to epoch 1998. In this 
Section, we model the dynamical evolution over the period 1912 - 1998, trying to find a model 
which gives the end of eclipses in 1950. 

In general, our model requires 25 input parameters, which we list here in two groups: 

QKa,, lAa, ^Aa, ^Aa! QAb, ^Ab; ^Ab, ^Ab (40) 

MAa, i?Aa, M Ah , R Ah ; Pa, eA; 

M B , Pab, eAB;«H, Pk, aj, Pj • (41) 

However, the A binary is sufficiently wide (P~ 14 d) that, provided its eccentricity (or more specif- 
ically its perihelion separation) does not vary, perhaps intermittently, by a substantial factor, tidal 
friction should be quite unimportant. More helpfully still, the Q-dependent distortion terms that 
determine X, Y, Z in equations (1) - (4) are unimportant compared with the third-body terms 
(components of the tensor S) in these equations, so that all of the quantities listed in (40) are 
negligible. The radii Ri and the angles defining the direction to the observer, are also 

unimportant for the orbital evolution, although they matter for the eclipses, and the date of their 
cessation. Although the Qi have little influence on the orbit, they do have a marked effect on the 
spins of the stars, because of the couple they cause - as we mention briefly below. 

This reduces our significant input file to the 13 quantities listed as (41). Of these, eA, Pa, cab, Pab 
are well or very well determined at epoch 1998 (TS00). Although eA may have (indeed will have) 
changed since epoch 1912, the other 3 quantities can be supposed constant. This is because, in our 
model, 

(a) at the level of the quadrupole approximation for the perturbing force of the third body, the AB 
orbit is exactly constant, and 

(b) the perturbing force on the A orbit, in the absence of tidal friction, is a potential force and 
hence does not supply energy to the A orbit when integrated over an approximately Keplerian 
orbit; this means that the semi-major axis cia, and the period Pa, are constant even though eA, h\ 
are not. 

Thus among these four quantities only e^(1912) must guessed - and ultimately solved for on the 
basis that in 86 years the 1998 value (0.136) must be reached (in conjunction with other constraints). 
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Similarly the 3 masses are constrained by but not uniquely determined by 3 observed mass- 
functions at epoch 1998. We need the two 1998 inclinations iab^a- TSOO estimated the latter 
on the basis that (i) za(1912) = 87.6° is known from the eclipse analysis, (ii) those data implied 
that i\ must have been 81.6° in 1950, when eclipses ceased, and (iii) ia has been decreasing at a 
constant rate since 1912. We find that in general the rate of change of ia is not very constant, and 
so we make a guess at the 1998 value of ia, which has of course to be consistent with the value 
that emerges from the calculation. The starting value aj is just aj = za(1912) = 87.6°. 

TABLE 2 

We also have to know or guess «ab- This however is constant in time since H - see (a) above - 
and J are constant vectors in space, even though their components in the e, q, h frame vary as the 
frame itself rotates. Dotting the vector H = (sin an cos fin, sin an sin /?h, cos an) into the vector 
J = (sin a j cos (3j , sin a j sin /3j , cos a j ) we have 

cos iab = cos qh cos qj + sin an sin aj cos(/3h — Pj) • (42) 

We therefore have to know or guess an and (3n — Aj in 1912, aj being known. 

TSOO's analysis of the ~ 1912 light curve gave an inclination of 87.6°, as mentioned above. 
They also obtained the longitude of periastron uj\ p , related to /3j by equation (18). The (3j from 
TSOO's 1912 light curve is 121.6° (see their Table 6, giving ui p , and differing slightly from their 
Table 5 value for reasons which they explain). However, this value is based on the assumption that 
eA is constant, and we find that generally it is not. The most significant orbital quantity which is 
given by the light-curve analysis, as TSOO explain, is the departure A<3?n = —0.072 of eclipse II 
from phase 0.5 relative to eclipse I. For moderate eccentricities, 

— A$n « eACOs^ip = -eAsin/3j = -0.1128 . (43) 

For eA = 0.136 this gives the value of /3j mentioned above. But in our best near-solutions we 
usually find eA increasing, i.e. it started in 1912 with a smaller value. Evidently it cannot have 
been smaller than 0.1128. Our preferred starting value is 0.115, and this implies /3j = 101.2°. A 
value above rather than below 90° is preferred, because TSOO's value of eAsina;i p = — eACOs/3j, 
though substantially less well determined, is fairly definitely positive. 

It may seem rather unsatisfactory that our preferred ca(1912) = 0.115 is very close to the 
minimum value 0.1128 inferred from eclipses. However, what can be seen as 'special' about the 
system is rather the fact that the 1998 value of /3j is extremely close to 90°: 91.7 ± 0.6 (TSOO, 
their Table 2, giving u\ p = 178.3°). Such a value, viewing the system almost exactly along the 
latus rectum, favours the maximum departure (A<&n) of the secondary eclipse from phase 0.5. If we 
imagine, going backwards in time, that eA does not change, then we are driven to postulate a rather 
large change in to the TSOO value 121.6°, to allow for the fact that the eclipses were substantially 
closer than this maximum value in 1912. What we conclude here is that less apsidal motion was 
necessary, because the eccentricity was a little smaller in 1912. We would quite generally expect 
that eccentricity changes on the same timescale as apsidal motion. Both timescales are dictated 
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primarily by the coefficient C in equation (15), if as for SS Lac only the third-body perturbation 
is significant. 

Our guess at the initial value of e& therefore provides us, from eclipse data, with a starting 
value for (3j via equation (43). We already have aj = i\ = 87.6° from TSOO's light curve data. 
We have to make two further guesses, at the angles ckr and (3u which in 1912 gave the orientation 
of H. 

To summarise, of the 13 quantities we need to start with in 1912, 4 are known directly from 
observation: these are Pa, Pab, &ab from the 1998 radial velocity curves and aj = «a from the 1912 
light curve. If we then guess the following 4 quantities - za in 1998, and 3 starting values eA, «h and 
Pb in 1912, we can work out the remaining 6 starting values from the following 6 observationally 
determined quantities: 3 mass functions from the 1998 radial velocity curve, and 2 fractional radii 
and the phase lag A<I>n from the 1912 light curve. Having integrated the equations for the 86 yr 
timespan, we then have 4 further pieces of observational data to constrain our 4 guesses. Three of 
these are eA and (3j = 270° — lu\ p in 1998, and the cessation of eclipses in 1950. We determine a 
theoretical Te, the time of cessation, as the average of the two times after t = (1912) at which 
the two series of eclipses (primary, and secondary) stopped. The observational value to match is 
Te = 38 yr. The fourth and last constraint on the four guesses is that the value for %a in 1998 
should equal the value guessed in the first place. Table 2 lists the values used, taken from TS00, 
and also lists our approximate solution. 

Table 2 groups parameters under 'observed', 'guessed', and 'computed'. All the observational 
data are taken from TS00. Our 'guess' was based on a preliminary eyeball search of parameter 
space, and refined by trial-and-error. 

Since the differential equations are non-linear, there is no guarantee either that a solution 
satisfying all the constraints exists, or that if it does it is unique. However our very brief search 
located one quite accurate solution with a fairly modest inclination between the orbits (29°), and 
a rather more extended search suggested that there were unlikely to be any other solutions, except 
possibly at high inclination where the behaviour can become rather chaotic. Another possibility, 
which we have not explored, is that the orbital inclination has decreased from 92.4° rather 87.6° in 
1912. 

FIG 5 

Figs 5a and 5b illustrate some aspects of the model. Fig 5a follows the orbital evolution for 
1912 - 1998, and Fig 5b follows it for slightly over 3000 yr. Date is plotted horizontally, and phase 
(from to 2, so that two complete cycles are shown) vertically. In Fig 5a, eclipses occurred within 
the long cigar-shaped areas labelled as I (*1 = Aa eclipsed by *2 = Ab) and II (Ab eclipsed by 
Aa). Eclipse I was slightly deeper and narrower, in 1912 (TS00). Phase in this figure is measured 
from periastron. Also shown as (i) and (ii) are the phases of points determined by equations (32) 
and (34). 

In Fig 5b, the same information is given for a much longer timespan: 1912 - 5250. Small leaf- 
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shaped patches now indicate the episodes of eclipses, and the curves indicating the phases (i) and 
(ii) slope generally downwards because of apsidal advance (and other rotation of the orbital frame). 
It can be seen that the next series of eclipses is not to be expected until about 2500. Epochs when 
eclipses take place appear to be separated alternately by a long interval and a shorter interval, and 
unfortunately we seem to be entering a long interval. Each period of eclipses lasts about a century. 

The fact that the phases of (i) and (ii) decrease (mostly, but not always) relative to the phase 
of periastron, as seen in Fig 5b, is of course due to the motion of the orbital frame. If this motion 
were just apsidal motion, i.e. if the major axis were rotating only about the angular momentum 
axis, we would have a relatively simple relation between the anomalistic period (periastron to 
periastron) and the sidereal period (successive passages through a plane fixed in an inertial frame 
and containing the CG), and lines (i) and (ii) in Fig 5b would have a constant slope. But with 
precession as well, the relation between the two periods can be rather complex. From the point 
of view of our simple model it is the anomalistic period that is 'basic': if the perturbing forces, 
however many of them, are conservative, the anomalistic period P as given by equation (16) is a 
constant at our level of approximation. 

The period of precession of the inner orbit is about 1000 yr, and the inclination to the line 
of sight (aj) oscillates between about 47° and 105°. The inclination of the inner orbit relative to 
the outer oscillates by only about 1°. It is inherent in our level of approximation that the inner 
angular momentum should be small compared with the outer, and unfortunately this is hardly true 
for SS Lac; but the fact that the inner orbit oscillates so little may nevertheless make the solution 
reasonably valid. 

We do not discuss the accuracy of the input data, and of our fit, in detail, for four reasons: 

(a) TS00 discuss fully the accuracy of the observational data. We have used only values which 
are independent of their assumptions that (i) ex is constant, and (ii) aj = i\ and /3j = 270° — uj\ p 
change at constant rates. All the standard errors are less than 1%, except for eAB (13%), A$n 
(6%), /b (3%), R/a (3%) and Te (3%); eAB only appears in C, equation (15), and in a very 
non-sensitive way. 

(b) We have zero degrees of freedom - 4 constraints to satisfy, with 4 unknowns - and so if we can 
find a solution at all it will be exact, to the extent that the data is. A hypothetical problem is that 
there might be functional dependences among the constraints, but the fact that our eyeball search 
converged very rapidly suggests there are not. Varying each of our 3 guessed angles by 1° usually 
gave a much worse fit, and so did varying eA(1912) by 0.001. Hence we believe that the guesses 
are right to about this level of accuracy. 

(c) By defining the computed value of Te as the average of the two times at which the two series 
of eclipses disappear, we make it a discontinuous (stepwise) function of the input, and cannot 
therefore differentiate it smoothly. We could develop a more sophisticated definition, but this 
seems unnecessary in view of the rather good solution found by trial and error. 



-24- 



(d) The main uncertainty is possible systematic error, such as the possibility that equations (1) - 
(4) are wrong. We had hoped to find more constraints than unknowns, and so test the theory more 
rigorously 

We can however make some predictions which are testable: for example the inclination 
should decrease to 70° in 2011, and to 65° in 2039. This should produce a measurable change in 
the mass- function. The eccentricity should be currently approaching its peak of ~ 0.138, and so 
may not change significantly for about 30 yr, but should drop to 0.132 by 2040. It should reach a 
minimum of 0.09 in 2160. 

FIG 6 

Fig 6 illustrates two possible behaviours of the spin The 3 components in the instantaneous 
orbital frame are shown as functions of time. The stars were started, arbitrarily, with Q. = uj. If 
the stars were perfect spheres they would simply maintain constant (vectorial) spin in an inertial 
frame - tidal friction being negligible in this system - and so oscillate sinusoidally in the frame 
of the precessing inner binary. But because they have quadrupole moments, due partly to their 
spin and partly to their gravitational effect on each other, there are couples on them. In Fig 6a 
we used our default value of Q = 0.028 (n = 3 polytrope), and in Fig 6b reduced this to 0.01. We 
have tried other values of Q, and do not see any very simple relation between the size of Q and 
the amplitudes or other characteristics of the oscillations. Considering that the orbit precesses on 
a cone of half-angle 29°, it seems surprising that the rotation axes of the component stars (we only 
plot *1 = Aa) can turn by more than 90° in the course of ~500yr. 

It may be questioned whether the approximation that we make in this paper, that a star 
rotates with a unique as if it were rigid, is sustainable in circumstances where Q, is changing 
in direction by a large amount in a few hundred years. Tidal friction is the agency that we rely 
on to achieve this: provided that the structure of a star is not strongly dependent on the velocity 
field within it, viscous dissipation should ensure that a non-uniformly rotating star evolves towards 
its minimum-energy state (for a given angular momentum) of uniform rotation. Although tidal 
friction between the two components of system A in SS Lac is probably negligible, tidal friction 
within either Aa or Ab is expected to operate on the timescale t v i sc of equation (9), i.e. decades. 
Thus it seems quite possible that the star is indeed kept fairly near a state of uniform rotation, 
despite major changes in the direction of its rotation axis. 

V907 Sco (B9.5V + B9.5V; 3.78d, e=0; Lacy, Helt & Vaz 1999) is another system in which 
eclipses come and go, even more dramatically than in SS Lac. It eclipsed in the intervals 1899 
- 1918 and 1963 - 1986, and not in 1918 - 1963, or after 1986. Lacy et al. (1999) detected the 
third body also from the motion of the CG of the short-period pair, with P out = 99.3 d, e ou t~0. 
Unfortunately an analysable light-curve for this system, during its eclipsing phase, does not exist, 
for reasons mentioned by Lacy et al. (1999), and thus we have less rather than more data with 
which to test our model. 
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5. Discussion 

Our formulation of the combined effect of five different perturbations on a Keplerian orbit, 
while very simple, appears to give physically believable results in a number of cases. It is however 
not easy to find observational data on stellar orbits which will seriously test the model. A major 
uncertainty is the viscous timescale of a star. One can question whether anything so simple as a 
unique viscous timescale can be adequate. However the timescale estimated from first principles in 
the Appendix seems surprisingly reasonable for the radio pulsar system 0045-7319. 

The model gives a determination of the orientation of the outer orbit in the triple system SS 
Lac, and though not over-constrained by data currently available may be challenged by data that 
should be available in a few decades. However in this system the quadrupole distortions of the 
stars are sufficiently insignificant that only the third-body terms are being tested here. 

A potentially significant statistical effect is predicted on the basis of the combination of tidal 
friction with Kozai cycles. We have argued that if the outer orbit in a triple is moderately highly 
inclined to the inner, then the inner orbit is likely to be shrunk to a limiting value of only 2 or 3 
days, supposing that it 'started' at a longer period. This limiting period will depend on the outer 
orbital period, and also on the rotation rate of the stars, etc. Roughly, it is dictated by the fact that 
for substantial Kozai cycles we need C^Z, from equations (12) - (15). This would give a longer 
limiting period for systems with a longer P out than f3 Per. On some models of triple star formation, 
inclinations larger than 60° are as likely as not, and so we might expect a significant deficit of 
orbits above some value in triples, relative to those in binaries. Tokovinin (p.c. 1998) has noted 
that in his multiple-star catalog (Tokovinin 1997) the distribution of periods among spectroscopic 
binaries that are in triples tends to drop off above 5 d, whereas among those that are not in triples 
it continues to rise. 

We have not yet been able to incorporate in the model a satisfactory approximation for what 
we believe is a very important further perturbation to binary orbits: the effect of ML and AML 
such as is likely to be experienced by cool stars with active dynamos in their outer convection 
zones. If a star is subject to spherically symmetric ML, the mass lost carries of orbital angular 
momentum as well as spin angular momentum: the amount of the latter may be enhanced if there 
is magnetic linkage between the star and the wind out to some substantial Alfven radius. However 
ML is unlikely to be very spherically symmetric. Also some of the wind may well be accreted by 
the companion star; and furthermore during the accretion process there is often found to be further 
ML, and presumably AML, in the form of bipolar jets from the inner portion of the accretion disc. 
A variety of possible models for the ML/ AML process can be thought of, but the physical process 
that they attempt to model may be too dependent on the details of how the gas actually travels 
from one star either to the other, or to infinity, to admit even at first order a simple yet credible 
formulation. We hope to attempt this in the future. 

An example where this may well be important is the young and active binary BY Dra (KlVe 
+ KlVe; 6d, e = 0.5; Vogt & Fekel 1979). One of the two components shows rotational modulation 
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with a period of 4 d. This is too slow for pseudo-synchronism, which at such high eccentricity implies 
a rotation period of 2 d. A possible answer is that the component is in a state of transient equilibrium 
between magnetic braking, which would tend to slow it down, and pseudo-synchronisation, which 
would tend to speed it up. It is by no means improbable that these two timescales are comparable 
in this system. 

Both the concepts of tidal distortion and of tidal friction will we hope be testable with some 
three-dimensional numerical modeling of stellar interiors that is currently being developed: the 
DJEHUTY project. This project aims at applying existing and well-tested 3-D hydrodynamic and 
thermodynamic (but non-self-gravitating) grid-based algorithms to the self-gravitating situation of 
stellar interiors, using massively parallel hardware. Although the resolution currently aimed for, of 
~ 10 s cells, would not be enough to resolve the surface layers of tidally-distorted stars, it may well 
be adequate to resolve the interiors and so determine whether dissipation in convective cores may 
be an effective agent of tidal friction, as we suggest here. 

This work was undertaken as part of the DJEHUTY project at LLNL. Work performed at 
LLNL is supported by the DOE under contract W7405-ENG-48. 
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The analysis that gives (a) the extra forces due to the five perturbative processes listed at 
the beginning of §1, and (b) their effect on the orthogonal triad e, q, h, is largely taken from 
EKH98, except for the third-body perturbation which was described in KEM98, and GR which is 
well-known. However EKH98 contains a mistake of a factor of 2 in the part of the gravitational 
potential that is due to the distortion of the stars by their mutual gravitational interaction. We are 
grateful to Dr. R. Mardling for pointing this out. We list here the equations of EKH98 that have 
to be changed: (36), (38), (75), (76) and (97). In each of these, the last term on the right, i.e. the 
term which does not involve Q, should be divided by two. This has no effect on the overall analysis 
in that paper. Equation (12) of the present paper, based on equation (97) of EKH98, contains the 
correction. 



This preprint was prepared with the AAS IATprjX macros v5.0. 



-28- 



Equations (10) and (11) of the present paper are obtained from equations (88) - (96) of EKH98. 
In EKH98 expressions were given for the rate of change of the Euler angles giving the orientation of 
the e, q, h frame relative to an inertial frame. In fact the X, Y, Z terms given here are what emerge 
more directly from the analysis; though not given explicitly in EKH98 they can be recovered from 
the formulae for rates of change of Euler angles given there. 

The the timescale ip for tidal friction that we use in this paper has been redefined to be twice 
the value that was used in EKH98. 

A novel result of EKH98 was a determination, exact to the order we work to here, of the velocity 
field in a rotating star which, in the frame that rotates with the star, suffers a time-dependent tidal 
perturbation due to the presence of the other star. Time dependence can arise because the orbit is 
elliptical, and/or not in corotation with the star. We summarise the result here. 

In the equilibrium-tide model, we approximate that the density (as well as pressure) is constant 
on equipotential surfaces of the instantaneous gravitational field of the companion. This means 
that 

P = p( r *) > r * = r + ra(r)P 2 (cos 0) , (Al) 

where a(r) is a dimensionless function that gives the ellipticity of an equipotential as a function of 
distance from the centre. The angle is measured from the direction of d(i), the separation of the 
two stellar centres. Radau's equation 



r z myr) \r r z 



gives a(r), apart from a multiplicative factor which gives a = a± at the surface r = R±, where 

M 2 Rl 1 

Q is related to the classical apsidal motion constant fo: 

Let 3 1 

F(r, 0) = r 2 P 2 (cos0) = -(k.r) 2 - -r 2 , (AS) 

where k = d/d. Since d is time- varying, both a and k depend on t, the former because a oc l/d 3 - 
equation (A3). Then with p as a function of r* only, and r* viewed as a function of r,t, we obtain 



dp dp dr* 
dt dr* dt dr. 



^SdaF + 3aG\ = ^ ^ (Idd \ 
ir* \dt r r J r ar* \ a dt J 



where 

1 dF d\i 
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F and G are obviously orthogonal harmonic functions of degree 2. 
Now consider the velocity field given by 



Then 



r ar \d ot 



and we can see that equation (Al) is satisfied to first order, provided that 

dp(3 



dr 



a dp 
cki dr 



i.e. 







' dp 
a— dr 
pot\ J Rl dr 



1 r 

POU Jr. 



(A8) 
(A9) 

(A10) 



The lower limit in the integral comes from the boundary condition that the outer surface (p = 0) is 
a surface which moves with the fluid, so that the velocity must be finite there despite the vanishing 
density. The function f3(r) is determined unambiguously by the structure of the star, via equation 
(A3) determining a(r), and is well-behaved ((3 — > 1) for polytropic (0<n<5) surfaces as p — > 0, 
despite the apparent singularity there. 



Using suffices, 

3ai at \ 
Vi = — (3{r)sijXj 



where 



- ddt WJ 



dkj dh 



(Aii) 



The rate-of-strain tensor is now seen to be 



_ dvi dvj 3qi / (3' \ 

tij = dx~ + dx- = ~2~ y 2 "^ + ~i SikXkX i + s jkXkXi\ J . (A12) 

We square this and average it over an equipotential (which at this level of approximation can be 
taken to be spherical), to obtain 



J = 9aj 4 (> + \rf30 + lr 2 /3' 2 ) 



Now, 



ldd\ 



6h77 +2 — = ^ 
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dd\ 



dd 



ddt) ' " \dt) d 2 
and so the rate of dissipation of mechanical energy is 
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The parameters w(r),l(r) are the mean velocity and mean free path of turbulent eddies. The /in- 
dependent weight factor in parentheses in equation (A15) is what we call 7(r), and its average over 
the turbulent convective region of the star is the 7 of equation (9). The factor in square brackets 
in equation (A15) leads to a functional form of the tidal friction force, as it depends on (variable) 
separation d, which is the same as the result usually obtained by arguing that the tidal bulge lags 
the line-of-centres by some small fixed amount. Averaged over a Keplerian orbit, it gives V and 
W - equations (5) and (6), and those parts of the terms X, Y, Z in equations (10) - (12) that arise 
from tidal friction. The details are given in EKH98. 

Although a common approximation for a(r) is a oc r 3 , and it is commonly argued from this 
that, in effect, f3 oc r 4 and 7 oc r 8 , none of these approximations is at all reliable. EKH98 integrated 
two polytropic models, and two MS stellar models. In the n = 3 polytrope, it was found (EKH98, 
Fig 1) that a and (5 decrease by a factor of about 10, from the surface right to the centre. At 
the surface (5 is unity, and 7 therefore somewhat larger (~4). In the central one-third by radius, 
roughly the region of convection in an upper MS star, 0.01^7^0.03. We therefore feel that an 
estimate of a mean 7 ~ 0.01 is reasonable, for MS stars which are typically slightly more centrally 
condensed than an n = 3 polytrope. 

The approximation a oc r 3 is appropriate to the outer layers of a star, where the density is low 
compared with the mean density: in this case equation (A2) gives m! = 0, i.e. p = 0. It is easy to 
integrate equation (A10) by parts, in this case, and the central value of (3 turns out to be just the 
ratio of mean density to central density. On the other hand, a =const. gives (3 = 1 throughout. 
The truth lies somewhere in between, with a ~ const, near the centre and a ~ r 3 in the outer layers. 
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Table 1: Proto- Algol binaries with different initial periods. 



Po(d) 


^cnd (d) 


T ciTC (yrs) 


io (deg) 


3 


2.8 


3xl0 7 


99.8 


5 


1.9 


lxlO 7 


97.3 


10 


1.7 


8xl0 5 


95.0 


15 


0.8 


9xl0 4 


94.0 



Table 2: System parameters for SS Lac. 



par am 


obs guess 


par am 


obs 


comp 


par am 


obs 


comp 


Pab 


679 d 


A$ n 


-0.072 




P-Aa/ a A 


.0741 






.159 


Aj 




101.2° 


RAb/a-A 


.0715 




Pa 


14.416 d 


lAB 




75.7° 


a A 




44.7 R Q 


«j = iA 


87.6° 


/Aa 


2.56 M© 




Ra& 




3.36 i? 


aj' = i' A 


73° 


/Ab 


2.49 M 




RAb 




3.20 R Q 


eA 


.115 




0.22 M 






.136 


.138 


«H 


29° 






2.93 M Q 




91.7° 


91.6° 




37° 


M Ab 
M B 




2.85 M Q 
.798 M 


T E 
aj' 


38 yr 


37.7 yr 
72.9° 



Note. — / Aa = M Aa sin 3 4;/ Ab = M Ab sin 3 4;/ B = M B sin4 B /(M A a + M Ab + M B ) 2/3 . 
Primed quantities refer to epoch 1998; all others, apart from Te, refer to epoch 1912, or to constants. 
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Fig 1 - The evolution of log(P/|P|) (plusses), dj (asterisks), (3j (circles) and V TO tsmi (crosses) 
with time in a binary like the SMC radio pulsar 0045-7319. Each quantity is divided by the 
observational value listed in the text. At about 195 yr, all four quantities are within about 20% of 
their observational values. At this point, (5q, 2 = 111°, having started with an arbitrary value of 0°. 
Such quantities as P, e, SI and a^ 2 have not changed significantly in this short time. 
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(a) Evolution of e, Oh/Om, P/PO (b) Omq vs Ome 
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Fig 2 - The evolution of orbit and spin in a binary like the SMC radio pulsar 0045-7319. (a) 
Eccentricity (plusses), cosine of the inclination of the stellar spin vector O2 to the instantaneous 
orbital plane vector h (asterisks), and period relative to initial period (circles). The B star was 
started with spin inclined at an 2 = 135° to the orbit (cos 0:0,2 = —0.71), reached inclination 90° 
after ~ 5 x 10 5 yr, and was almost completely parallelised by ~ 1.7 x 10 6 yr. (b) The two components 
of B-star spin in the orbital plane, plotted against each other. The spin axis started at the top, 
left-of-centre, turned through ~1.25 rotations anticlockwise, then (at the time when the spin was 
exactly perpendicular to the orbit) reversed its motion to clockwise while spiralling in towards the 
centre. The evolution was speeded up by artificially decreasing the viscous timescale by a factor 
of ~ 500, to prevent the spiral being very tightly wound. The 'real' timescale would have required 
about 600 turns before reversal, (c) Timescales e/|e| (asterisks) and P/\P\ (plusses); also the 
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timescales \/Z (circles) and (crosses) of apsidal motion (all logged). The first two timescales 

are artificially shortened by ~500, as in (b). The last two timescales tend towards equality as the 
spin becomes parallelised, but precession due to non-parallel spin causes one to oscillate about the 
other, (d) The precessional timescales 1/V X 2 + Y 2 (plusses), and the timescale of rate of 

change of inclination to the line of sight (both logged). The many cusps in the latter are due to 
the fact that the inclination was oscillating between two values. 
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Fig 3 -The Darwin (D) and eccentricity (E) instabilities. Eccentricity (plusses), orbital fre- 
quency uj relative to its initial value (circles), the degree of asynchronism, log(f2/cj) (asterisks), and 
the ratio of spin to orbital angular momentum, \og{IVL / fj,h) (crosses). *1 is a neutron star, and 
*2 a massive, partly-evolved MS star. The orbit has P = 6d, e = 0.1 to start with, (a) Initially 
0,/uj = 2. (b) Initially = 0.7. In (a) the system starts both D-unstable and E-unstable. 
Eccentricity and asynchronism grow, but the periastron separation remains large enough to avoid 
collision. Once the orbit has widened it becomes stable to both processes, and settles down. How- 
ever nuclear evolution (neglected) would cause problems before 10 7 yr. In (b) the orbit is E-stable 
and slightly D-stable to start with. But as the orbit and star gradually spin up the orbit's angular 
momentum goes down, while the star's goes up, leading to D-instability in about 5000 yr. After 
that asynchronism increases, and the stars collide in about 19000 yr. 
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(a) short term (b) medium term (c) long term 




Fig 4 - Evolution of eccentricity (dots) and log P (thick line) in the inner binary of a 'proto- 
Algol' triple system. The initial orbital parameters are ((2.5 + 2 M@; 5 d, e = 0.1) + 1.7M ; 679 d, 
e = 0.23; «h = 97.5°). (a) The first 2000 yr, showing somewhat truncated Kozai cycles; (b) the 
first 10 6 yr, showing the orbit settling towards a nearly constant but slowly decreasing eccentricity; 
(c) the first 10 7 yr. By 10 7 yr, e~0, P~2.1d, and «h = 100°. Some apparent structure in the 
eccentricity variation in (b) is due to beating between the data-plotting frequency and Kozai-cycle 
frequency, which can be commensurable. 
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(a) phases of eclipses 1912 — 1998 (b) phases of eclipses 1912 — 5250 




Fig 5 - (a) The phases of eclipses of SS Lac as computed here for 1912 - 1998. Two whole cycles 
are shown on the vertical axis. Phase zero is periastron. Eclipses occur in the narrow cigar-shaped 
areas centred at phases 0.24 (eclipse II, *1 in front) and 0.81 (eclipse I, *1 behind), and ending 
at about 1950. Also shown are two phases labelled (i) and (ii). The first is where *1 crosses the 
plane containing the line of sight and the orbital axis, *1 in front - equation (32) - and the second 
where the radial velocity of *1 relative to the CG of the inner binary is zero and decreasing (i.e. *1 
behind) - equation (34). The slight effect on phase of the orbital motion of the inner binary within 
the outer binary has been ignored, (b) Same as (a) but for the interval 1912 - 5250. Regions 
of eclipses are now leaf-shaped. The sloping lines can be identified by comparing the left-hand 
edge with the whole of (a). The slopes indicate that periastron is, on the whole, advancing, but 
occasionally retreats because of precession. 
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(a) Components of spin: Q=0.028 (b) components of spin: Q = 0.01 




Fig 6 - The three components of the angular velocity of *1 in SS Lac as functions of time: 0,^ 
(circles), Oi e (plusses) and f2i q (asterisks), (a) Qi = Q\ = 0.028; (b) Q\ = Qi = 0.01. In both 
cases, the system was started with parallel corotation. 



